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Complex mesoscopic systems play increasingly important roles in modern science - from under¬ 
standing biological functions at the molecular level, to designing solid-state information processing 
devices. The operation of these systems typically depends on their energetic structure, yet probing 
their energy-landscape can be extremely challenging; they have many degrees of freedom, which 
may be hard to isolate and measure independently. Here we show that a qubit (a two-level quan¬ 
tum system) with a biased energy-splitting can directly probe the spectral properties of a complex 
system, without knowledge of how they couple. Our work is based on the completely-positive and 
trace-preserving map formalism, which treats any unknown dynamics as a ‘black-box’ process. This 
black box contains information about the system with which the probe interacts, which we access 
by measuring the survival probability of the initial state of the probe as function of the energy¬ 
splitting and the process time. Fourier transforming the results yields the energy spectrum of the 
complex system. Without making assumptions about the strength or form of its coupling, our probe 
could determine aspects of a complex molecule’s energy landscape as well as, in many cases, test for 
coherent superposition of its energy eigenstates. 

I. INTRODUCTION 

Measuring the properties of complex systems at the 
mesoscopic level is extremely difficult and of great im¬ 
portance. By definition, these systems are too small to 
admit an effective statistical description, but too large 
for their dynamically relevant degrees of freedom to be 
straightforwardly isolated. This is especially true at the 
classical-quantum boundary, where quantum effects are 
still important, but understanding their dynamical and 
energetic consequences can become highly non-trivial due 
to the size of the system [1]. Such systems play a vital 
role in the function of living organisms [2-8] , and increas¬ 
ingly form the basis for quantum technologies [9]. Over 
the last few decades, much effort has gone into developing 
experimental techniques to characterise and quantify the 
salient features of complex systems in a variety of con¬ 
texts. For example, pump-probe and multidimensional 
spectroscopy are used to isolate electronic contributions 
to the energetic structure of complex molecules [10, 11[ 
and solid state systems [12[, while techniques utilising 
Raman scattering are used to determine vibrational spec¬ 
tra [13, 14[. In the same spirit, there have been several 
recent attempts to quantify quantum coherence in com¬ 
plex systems [15-19[. 

To learn anything about a system, it must first in¬ 
teract with a probe, which induces open dynamics for 
both the probe and the system. From the perspective 
of the probe, the complex system acts as a black box, 
which evolves the probe’s quantum state in a generally 
non-unitary way. Information about the complex sys¬ 
tem must then be determined by studying this black-box 
dynamics. The mathematical tools that describe open 
quantum dynamics operationally have a rich history [20[ . 

Foremost among them is the completely-positive trace¬ 


preserving (CPTP) dynamical map [21, 22[, which de¬ 
scribes any transformation between two temporal points 
in full generality - the map At-:o relates the probe’s den¬ 
sity operator at the end of the evolution, Pt, to that at 
the beginning, pQ. The map contains all the dynamical 
information that is experimentally accessible by making 
measurements on the probe only, and we can ask what 
one can infer from it about the Hamiltonian of the com¬ 
plex system Hs and its interaction with the probe V. 
Stinespring’s dilation theorem tells us that any CPTP 
map stems from unitary evolution of the probe and sys¬ 
tem together [20[: 

Pr = A^.o[po] = trs [e"*^*°‘ ""po ® Ps '^] , (1) 

where i?tot = H-p -\- Hs + V, Hp can be identified as the 
Hamiltonian of the probe, ps is the initial state of the 
system (which we will always assume to be uncorrelated 
from the probe initially), r is the evolution time and we 
have set h = 1. However, this relationship is not unique; 
in fact, without additional assumptions about the nature 
of the system, there is little that can be directly deduced 
from At-:o - the dynamics of the probe alone - about ps 
and Vps = Hs + V. 

Here, we develop a new method for probing, in prin¬ 
ciple, any uncharacterised complex system, be it a large 
organic molecule or some exotic meta-material, using the 
simplest quantum probe - a qubit, or two-level quantum 
system. Several schemes already exist which use qubits 
to probe larger systems: impurities in Bose-Einstein con¬ 
densates can be used to determine phase [23[, phonon 
spectra [24[ and temperature [25[; in cavity and cir¬ 
cuit quantum electrodynamics, effective two-level sys¬ 
tems have been proposed as probes for environmental 
noise spectra [26, 27[ and witnesses for strong coupling 
[28[; the use of single spins to precisely map complex 
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electric and magnetic fields [29] and environmental co¬ 
herence [30] has also been proposed. However, in each of 
these cases, the analysis relies on a particular form for 
the probe-system Hamiltonian Vps- Our result differs, 
and is more generally applicable, in that it applies even 
when Vps is completely unknown, though in many sit¬ 
uations it may be limited by an inability to control the 
probe sufficiently. 

It should be noted that Vps necessarily includes those 
degrees of freedom which would usually be considered as 
an environment to the system of interest. As long as 
these do not strongly couple directly to the probe, their 
effect will generally be to smear out the energy levels of 
5. 

The scenario we consider in this article is illustrated 
in Fig. 1. We allow a classical control field with variable 
magnitude (denoted by A) to modulate the probe’s en¬ 
ergy splitting in a basis of our choosing, parameterised 
by the angles 9 and (j), leading to a Hamiltonian for the 
probe 

H-p = (2) 

where cr^e,^) = sin 6* cos ^ (Ja, -I- sin6*sin^CTy -I- cosflcr^ is 
the Pauli operator that determines the basis in which 
the bias is applied. Moreover, we are free to choose the 
input state of the probe and, after a time r, to make a 
measurement of the probe state. 

Our setup is inspired by [18], where Markovian master 
equations are studied with an added control field. Within 
this minimal scenario we show that it is possible to infer 
information about the unknown system, including prop¬ 
erties of its state and spectrum, directly from the statis¬ 
tics of the final measurement of the probe. We find three 
universal regimes ordered in powers of A, each of which 
witnesses different properties of the system. In particu¬ 
lar, we show how to measure time-correlation functions 
for the system operators, which indicate the timescales 
over which the system behaves quantum mechanically. 
The simple probe that we describe here could be de¬ 
ployed to investigate a large variety of systems, and the 
more that can be assumed about a particular system, 
the more specific the information than can be obtained 
through our method, as we demonstrate with two exam¬ 
ples - one involving spin systems, the other concerning 
vibrational degrees of freedom. 

II. PERTURBATIVE EXPANSION IN THE 
PROBE ENERGY SPLITTING 

Let us denote the prepared state of the probe as |a) 
(our results can be straightforwardly extended to mixed 
states by linearity) and define the basis in which we 
measure by the state |/3). The experimentally accessi¬ 
ble quantity (from which the dynamical map A 7 .:o can be 
recovered) is the transition probability 

Pp-.a = (/3|AT-:o(|a)(a|)|/3) • (3) 
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FIG. 1. (Color online) The controlled two-level probe 
scenario: The properties of an uncharacterised complex sys¬ 
tem 5 can be inferred from the way in which it influences 
dynamics of a two-level probe V, with which it interacts. By 
varying a controllable bias of magnitude A, more information 
about the system becomes accessible, including aspects of its 
state and spectrum. 

By determining this quantity for different preparations 
and measurements as the control parameters are varied 
- as depicted in Fig. 2 - information about the system can 
be recovered. Furthermore, as we will show, much of this 
information can be recovered using preparations, mea¬ 
surements, and the control field in a single basis, mean¬ 
ing that our scheme could be used even with probes that 
have a preferred measurement basis. 

The method presented here is general; however, its 
practicality will be determined by our ability to increase 
the energy splitting of the probe, as well as to make 
preparations and measurements. In practice, this may 
require that the probe is prepared independently (away 
from the system) and brought into and out of contact 
with the system in a repeatable manner, such that the 
two are not interacting when preparations and measure¬ 
ments are made. Alternatively, the probe can be placed 
next to the system and prepared via control operations 
on the probe alone, which must act on a timescale faster 
than that set by the interaction with the system. 

Combining Eqs. (1) & (3), the probability of finding 
the probe in the state |/3), after the open dynamics takes 
place, can be written as the expectation value of the pro¬ 
jector \/3){(3\ 0 1 with respect to the time evolved density 
operator: 

Pfs-.a = tr [(|/?)(^| 0 l)C/(|a)(a| (g) psW^ , (4) 

where U = exp {—iHtotT} is the time evolution op- 
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Davies in order to derive a weak coupling master equa¬ 
tion [31]. 

To simplify the subsequent calculation, we can always 
decompose V-ps in the basis defined by the eigenstates of 
the control Hamiltonian, {Itto) , Itti)}, 

Vvs =ko)(7ro| 0 Ag -k IttoXttiI (g) 

+ ki)(7ro|(g)H=i-k|7ri)(7ri|(g)Hi, (5) 


FIG. 2. (Color online) Probe evolution as a black box: 
(a) A two-level probe is prepared in state |a); (b) the probe 
is acted on by CPTP map At:o(A, 6, <(>), which has a set of 
controllable parameters corresponding to the magnitude and 
basis of an imposed splitting of the probe’s energy levels H-p = 
(c) the probe is measured in the basis defined by 

state I/?). 


erator for the joint probe-system. In general, due 
to the unknown nature of Vpsi it is not possi¬ 
ble to directly calculate U. However, by rewrit¬ 
ing it as C/ = exp{-i (iAcr(e_,^) (g) 1-I-Pps) r} = 
exp{—i (A 0 1 + V-ps/^ 'At}, one can see that, in 
the limit of large A, Vps acts as a perturbation to the 
measurement Hamiltonian, albeit evolved to an effective 
time A r. We can then move into the interaction picture 
and and write down the transition probabilities as an ex¬ 
pansion in 1/A. This relative scaling of the probe and 
interaction Hamiltonians is related to that proposed by 


with Aq = Aq^ and A^ = A^K This has the advantage 
of further specifying how different system operators cou¬ 
ple to the probe. By varying the angles 9 and (/, and 
hence the basis in which the bias is applied, different sys¬ 
tem operators can be selected to couple to different probe 
states (see Appendix A for details). From this point on 
we drop the explicit basis dependence for the system op¬ 
erators {e.g. Ag = Ag). 

If we move into the interaction picture with respect to 
Ho = 0 1 + (1/A)|7rg)(7rg| 0 Ag -b (1/A)|7ri)(7ri| (g 

Ai, then the corresponding time evolution operator has 
a perturbative expansion. The resulting Hamiltonian is 
Ht = exp{it}|7rg)(7ri|gi3t-|-exp{-ft}|7ri)(7rg|gi3|, where 
Bt = exp{iAgt/A}H exp{—fAit/Aj. 

In the interaction picture, the time evolution opera¬ 
tor can be expanded in a Dyson series as t/j = 1 + 
(1/*A) f ds Bs-(1/X^) Jo' ds /; ds' BsH,, +0 (l/A^) , 
which in turn leads to an evolution equation for the den¬ 
sity operator pt = Ut po U} ^ where 


rXt 


rXt 


Pt =Po + ^ / ds[Ffs, Po] - ^ 


ds 


ds^ ( HsHs'Po + pqHs'Hs ) + 


rXt 


pXt 


ds / ds' HsPoHs' + 


( 6 ) 


The Schrodinger picture density operator is given by pt = exp{—lATFgt} pt exp{iAi?gt}, and the transition probability 
after evolving for time r can be expressed, up to second order in 1/A, in terms of interaction picture quantities as 
(see Appendices B & C for details and for the the conditions of validity of the perturbative series) 


Pp-.a =tr 


1 

V 


1 0 

bobl 


+ alaxhoh\C,^°'>{X) + j oou} + bobi (|aoPCo^^(A) - |aip^J^^(A) 


(|aoPer(A) - |ai|2ef)(A)) +aSai6o6}C^"HA) 


+ c.c. + 0( - 


(7) 


( 8 ) 


where qp,^ = j&oaSr + l^ia}| , Ofc = (a| tt^) , bk = {/3\ TTk), and 

^(0) (y) ^ g*Ar ^ ^(1) ^ (^^^Aor^ ^ , 

^(1)(A) = - e'^^ B^)^ , ^^(A) = - e*^^ , (9) 

= {BB"!)^ + - e*^^ BB'f)^ , (10) 


^+XA) = -b - e*^^ {Bh 


^Be- 


%■ 


( 11 ) 


Above {■)g =tr[-p 5 ] and the functions {^(®)(A),^^*^(A)} all depend on r, Vps and ps', their magnitudes remain 
















4 


constant with A. written out explicitly 

here, as it never appears at leading order in A. It should 
be noted that Eq. (7) is linear in the choice of prepara¬ 
tion and measurement state, and can thus be straight¬ 
forwardly generalised to mixed initial states pa and gen¬ 
eral measurement, positive operator valued measure, ele¬ 
ments by making the replacements Uk = tr[pa|7rfe)(7rfc|] 
and bk = tr[n^|7rfc)(7r/j|] respectively. These more general 
choices would introduce uncertainties to the procedure 
outlined below. 

In Appendix B, we discuss in detail the conditions 
required of A for Eqs. (7-11) to accurately describe 
the transition probabilities. In short, we require that 
^ I (iil^l^o) I) where |ji) and |fco) are eigenstates of 
Aq and Ai respectively that have support in - or are eas¬ 
ily reached from - the initial state ps . It is also necessary 
that the probe is not resonant with any energy splittings 
in the system that can be reached within three appli¬ 
cations of the perturbation Hamiltonian Hs (or at least 
that such resonances satisfy | {ji\ B |fco) | l/(''^^'^))- 

There is an important further restriction that depends 
on the leading order A dependence in Eq. (7) and which 
generally relates the value of A to the evolution time r. 
For zeroth and second order dependence, the constraints 
are 

^ » E (jol B^B\k^){k^\ps Ijo) I (12) 

jk 

and 

JK 

respectively, where 

' \r E‘=El ' ' 

Here and Ej are the eigenenergies of Aq and 
Ai respectively., i.e., Aq = Y.j Ej\jo){jo\ and Ai = 
The same expressions should also hold 
for 0 —>■ 1 and B B^ . When the transition probability 
is proportional to 1/A to leading order, there is a similar 
constraint that is independent of A: 

tr{Bps} > (jol B'fB\ko){ko\ps |jo) I, (15) 

jk 

and the same with 0 —?► 1 and B ^ B^. 

With the exception of the latter, first order case, 
and practical limitations notwithstanding, sufficiently in¬ 
creasing the value of A will always lead to probabilities 
consistent with Eq. (7). Even though the signature oscil¬ 
lations of pp:a with A arise for energy splittings inconsis¬ 
tent with the constraints in Eqs. (12) & (13), their valid¬ 
ity could in principle be checked by further increasing A 
and again measuring the transition probabilities. Once 
the amplitudes of oscillation have converged, one can be 
sure that Eq. (7) correctly describes the dynamics. 


III. INFERRING PROPERTIES OF THE 
SYSTEM 

We can see from Eq. (7) that the dominant term at 
large A depends on the choice of |a) and \/3): in gen¬ 
eral, zeroth order terms dominate, but when either the 
preparation or measurement is made in the {Itto) , Itti)} 
basis, the leading order contribution is proportional to 
1/A. When both |a) and \(3) are eigenstates of H-p, then 
only terms that are second order in 1/A survive. In each 
case, however, the probability for large A - subject to the 
conditions discussed in the previous section - can always 
be written in the form pi},a — Qp-.a + 2Re{A}, where the 
dependence of X on the choice of preparation and mea¬ 
surement is summarised in Table I. The real part of X is 
an oscillatory function of A: 

2Re{A} ~ 77 -I-£)cos(At-I-(/ j), (16) 

with frequency r, an amplitude D, phase ip, and a shift 
rj. All of D, ip and 77 depend on the system state and 
Hamiltonian. 

This method for inferring properties of the system 
bears analogies with scattering theory. In both cases, 
a known wave function is used to probe the structure of 
an unknown object, resulting in phase shifts and ampli¬ 
tude/probability changes. By contrast, in the case inves¬ 
tigated here, the probe is always coupled to its environ¬ 
ment, so one has to find a way to switch off the interac¬ 
tion or, alternatively, to uncouple the two systems. Also, 
rather than looking at typical scattering quantities, such 
as the cross section, we focus on other structural proper¬ 
ties of the system and the interaction, such as its survival 
probability in the initially prepared state. 

In Fig. 3, we plot pp,a. as a function of A for several 
different choices of preparation and measurement (with 
a randomly chosen Vps and ps for a four-dimensional 
system). The probability can be seen to converge on a 
single oscillatory component with an appropriate enve¬ 
lope function. By characterising these oscillations, the 
values of the expectation values appearing in each of the 
terms of Eqs. (8-11) can be determined. If we further 
allow to vary the length of time for which the probe in¬ 
teracts with the system, then we can determine quantities 
such as and B ^as a func¬ 

tion of r. We now discuss which elements of the system’s 


X 

k)^{ko),ki)} 

k)e{ko),ki)} 

|a) i {ko) > ki)} 
|a) = kfe) 

flSai&ok k°kA) 
(-l)"bS 6 iCi''(A)/Ae 

±ooaK'^k^)/A 
(-l)"+kf (A)/A^ 


TABLE I. Leading order A dependence of pp-.a- The 
functional dependence is always oscillatory, with amplitude 
and phase that depend on the quantity X given in Eq. (16), 
and an envelope function which depends on the order. The 
choice of preparation (|a)) and measurement (|/?)) determines 
X, and hence the dependence of pp-,a on Vvs and ps- 
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FIG. 3. (Color online) Transition probability. Typical be¬ 
haviour of pp-a as A is varied from small to large (for ref¬ 
erence the matrix elements of Vvs are sampled such that 
I (io(i)l B |fco(i)) I ^ V^). The results shown are for a system 
comprising two qubits; Vvs and ps were sampled randomly, 
such that they were effectively unknown prior to the simu¬ 
lations. Different panels show effects at different orders for 
large A, corresponding to different choices of preparation and 
measurement. Panels (b) and (c) both correspond to first or¬ 
der effects, but with different choices of preparation. Solid 
red represent the values of pp-a, and bounding envelopes are 
shown with dashed blue lines. 


energy spectrum and state can be recovered from these 
measurements, for more detail, see Appendix D. 

When neither the preparation nor the measurement is 
in the basis of H-p, the terms of Eq. (16) become rj = 0 
and He"^ oc with a coefficient that de¬ 

pends on |a) and |/3). Thus, the zeroth order term mea¬ 
sures how differently the two eigenstates of the probe 
couple to the system. Oscillations that depend on system 
operators Aq and Ai can be attributed to ‘Lamb-shift’ 
terms in the probe’s Hamiltonian: an effective level split¬ 
ting arising from its interaction with the system. The ze¬ 
roth order becomes trivial, i.e., the oscillations no longer 
depend on the system when Aq = Ai, i.e., when there 
is no interaction term proportional to Hp (in this case 
Aq = Hs, the Hamiltonian of the system). 

When either |a) or |/3) (but not both) are eigenstates 
of Hp and Eq. (15) is satisfied, all zeroth order A de¬ 
pendences vanish and terms proportional to 1/A domi¬ 
nate the large-A behaviour. From Table I, it can be seen 
that there are three distinct cases in which first order 
effects dominate, each measures a different property of 
the system. In particular, the Fourier components of the 
quantities have frequencies equal to 

energy differences between eigenstates of Ak which are 
initially in coherent superposition. Furthermore, when 
the central value rj of A-oscillations dictated by or 
is independent of r for some choice of 6 and (p, then 
the system state ps is diagonal in the basis of the corre¬ 
sponding operator, Aq or Ai respectively. Remarkably, 
when Aq = Ai, oscillations at this order can therefore be 


used to test whether a system is in equilibrium. 

Finally, when both preparations and measurements are 
made in the basis of Hp, only second-order and higher 
terms survive in Eq. (7). In this case, the system quan¬ 
tities that can be inferred are all two-time correlation 
functions for the operators Bt- These two-time correla¬ 
tion functions (also known as bath or reservoir correlation 
functions) dictate how the system influences the probe’s 
dynamics in the absence of the control bias, quantifying 
memory effects and dynamical coherence [1, 20]. The 
second order effect can be related to the quantum Zeno 
effect [32, 33] as shown in [18]. The Fourier components 
of the correlation functions correspond to energy differ¬ 
ences between all eigenstates of Aq and Ai that have 
support in ps and those to which they couple through 
B. In other words, we get insight into the spectrum of 
Vps - for a closed system this will be a series of dis¬ 
crete levels, whereas a system which is in turn coupled 
to a large reservoir will have a continuous distribution. 
Lastly, the decay of A-dependent oscillations at the sec¬ 
ond order (and hence the correlation function) above a 
certain evolution time r signify a point beyond which sys¬ 
tem memory effects do not persist. Thus we can quanti¬ 
tatively measure the forgetfulness of a system, which is a 
generic property for large numbers of degrees of freedom. 

Qualitatively, this behaviour is generic for quantum 
systems - even when the system is further coupled to 
some environment. In this case the B operators will cou¬ 
ple subspaces of the larger closed system composed of S 
and its environment and, in many cases, there will be a 
continuum of energy spacings in the Aq and Ai opera¬ 
tors. However, there are several exceptional cases that it 
is worth mentioning in more detail: 

If the system S is classical and time-independent, 
H (X Ag oc Ai and the probe undergoes Hamiltonian 
dynamics; it must be considered as an isolated quantum 
system under the action of some potential. Since the po¬ 
tential is “classical”, it will exhibit no recoil effects, and 
as a consequence, the quantum dynamics of V will not 
be able to unveil any interesting quantum-like feature of 
(classical) S. 

The same considerations would apply for a “slowly” 
changing classical system S. In such a case the Hamil¬ 
tonian of the probe V becomes time-dependent, and our 
results would still hold with all the A’s and B’s replaced 
by scalar functions of time. 

On the other hand, a classical stochastic system S 
would induce (unital) CPTP dynamics for the probe, by 
virtue of the Stinespring dilation theorem: the dynamics 
would be governed by a Gorini-Kossakowski-Sudarshan- 
Lindblad (self-dual) master equation, with potentially 
time-dependent rates. Our method would then yield 
the two-time correlation functions for B, which is now a 
stochastic variable. For more general semi-group dynam¬ 
ics, there is always a consistent time-independent Hamil¬ 
tonian for a larger system that would produce the same 
behaviour. 
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IV. EXAMPLE SYSTEMS 

Of course, it is rare that a system is completely un¬ 
characterised. More often than not, there is an assumed 
Hamiltonian for the system, which may include unknown 
parameters. It is also common to assume that the sys¬ 
tem is in thermal equilibrium with its environment. In 
such cases, assumptions about the system can be used 
to extract more specific information from the correlation 
functions discussed above. We now briefly present two 
generic and commonly used models for complex systems 
as examples, namely those comprising collections of spins 
and collections of harmonic oscillators (representing the 
vibrational modes of a molecule). The former tend to be 
used when the relevant degrees of freedom are highly lo¬ 
calised [34], whereas the latter are a good approximation 
to systems with delocalised degrees of freedom close to 
equilibrium [35]. 


Spin systems 


As our first example, we will consider a system com¬ 
prising a small collection of spins - for example molecular 
nuclear spins or isotopic impurities in otherwise spinless 
solids (for example silicon [36]). The energetic structure 
of such systems is routinely studied using nuclear mag¬ 
netic resonance (NMR) spectroscopy, and characteristic 
spectra are used to identify different molecular species 
[37]. However, obtaining NMR spectra usually requires 
a large ensemble of sample systems and strong mag¬ 
netic field gradients. We will now demonstrate that our 
method could in principle be used to obtain spectra for 
single molecules using, for example, a tuneable magnetic 
dipole as a probe, without precise knowledge of the rela¬ 
tive location of the sample with respect to the probe. It 
should be noted that our approach is not unique in this 
respect; other proposals have been put forward for single 
molecule NMR spectroscopy [38] . A typical example of 
such a tuneable magnetic dipole would be a supercon¬ 
ducting qubit, such as a flux qubit, although such probes 
are often fairly strongly coupled to wider environments. 
In practice, this coupling leads to rather short coherence 
times for flux qubits - state-of-the-art Ti and T 2 times 
are ^ 40^s and ^ 80fj,s respectively [39] - limiting the 
energy resolution of reconstructed spectra. What follows 
should be viewed, in the context of current technology, as 
a proof of principle for our scheme, rather than a concrete 
experimental simulation. 

Spin systems in an external magnetic held Bq have a 
Hamiltonian comprising a Zeeman splitting term and a 
magnetic dipole-dipole interaction: 


.^spin 


~ • Sfc -1- 




St ■ Sfc, - 


47r|rfcfc/|3 

(Sk ' ^kk''){Sk' * I'/cfcO 
Irfcfe'P 


, (17) 



0.0 0.2 0.0 0.2 0.4 

Energy (neV) 


FIG. 4. (Color online) Simulated spin spectra: (a) A su¬ 
perconducting flux qubit probe V (with magnetic moment 
0.3 meV T”'^) is placed at four different randomly chosen posi¬ 
tions (i-iv) within a radius of 100 nm around a spin system S, 
which consists of four spins (with magnetic moments ~ /rjv) 
randomly arranged within a 0.01 nm sphere; a |B| = ImT 
magnetic field is applied perpendicular to the magnetic mo¬ 
ment of V. Probe positions are shown in the plane defined by 
B and the vector joining V and S. (b) Reconstructed spec¬ 
tra for the different probe positions. Red peaks indicate the 
spectrum that would be obtained with perfect measurement 
of the transition probabilities pi.,o for all values of A and r. 
The other curves are simulated reconstructions, where the ex¬ 
periment was performed only a finite number (~ 10®) of times 
for each choice of parameters, r was sampled evenly in 100ns 
steps for different total evolution times: 80 ps (dotted green); 
160/is (dashed purple); 2ms (solid blue). 100 different values 
of A were tested in each case. 


where rtk' = — r^/ and Sk, and fXk are the spin 

operator, position and magnetic moment of the fcth spin 
{fiQ is the permeability of free space). It is the differences 
between the eigenvalues of Eq. (17) - usually broadened 
by the snrrounding environment - that determine the 
positions of peaks in NMR spectra. For nuclear spins in 
small molecnles, the magnetic moment is on the order of 
the nuclear magneton /ijv and the nuclei are often only 
a few angstroms apart [37] . In the absence of an exter¬ 
nal held, this leads to typical energy splittings of aronnd 
lOpeV. 

The qubit probe itself is a magnetic dipole with mag¬ 
nitude /i-p, which is a function of the persistent cnrrent 
running through it; the orientation of this dipole mo¬ 
ment depends on the quantnm state of the qnbit. If we 
nse the two states of definite magnetic moment to define 
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the eigenstates of the probe’s ax operator, then we can 
write its interaction with the spins as 


H- 


Pspin 


=E 


flOfJ-VfJ-k 

47r|r-pfc|3' 


X [n-p -Sk- 


(18) 


{fip ■ r-pfc)(Sfe • r-pfe) 


where n-p is a unit vector along the probe’s magnetic 
moment when it is in the +1 ax eigenstate and is 
the vector joining the probe and the fcth spin. It is also 
possible to tune the energy splitting of the flux qubit 
in Ref. [40] between the eigenstates of cr^ such that the 
persistent current (and hence magnetic moment) does 
not change appreciably, that is it satisfies the conditions 
required of a probe for our scheme. 

Since only transitions between orthogonal probe states 
are of interest in this scenario, preparation of the probe 
can be performed using the same kind of projective mea¬ 
surement used to determine which of the two eigenstates 
of az the probe is in. For flux qubits, such measurements 
can be performed on timescales of a few ns - much shorter 
than the interaction timescale between probe and sys¬ 
tem - by biasing the energy of the definite current states 
(eigenstates of ax) and coupling it to, e.g., a Josephson 
bifurcation amplifier [41, 42[. 

In Fig. 4, we show a simulated reconstruction of the 
spectrum of a randomly positioned collection of spins in 
an external magnetic held, described by the Hamiltonian 
in Eq. (17). The plots in the Figure are the Fourier trans¬ 
form of the two-time correlation functions in Eq. (10), in¬ 
ferred using our scheme from measurements of the transi¬ 
tion probability between orthogonal states of the probe; 
the different curves correspond to different total evolu¬ 
tion times. Despite the fact that the probe is placed at 
different distances from the system, the same spectra are 
reconstructed, with greater or lesser accuracy. Only in 
those cases where the coupling is very strong, and the 
value of A necessarily very large, does the reconstruc¬ 
tion fail; this is due to the difficulty of measuring small 
probabilities in a finite number of experiments. As can 
be seen from the Figure, the energy resolution of the re¬ 
construction is limited by the evolution time of the joint 
probe-system; this is due to the comparatively long time 
scales associated with the spin system. In performing 
these simulations, we have assumed that other environ¬ 
ments do not play a significant role in the dynamics on 
these time scales. While this is a poor assumption for 
flux qubits with current technology, where 1 // noise lim¬ 
its coherence times as discussed above, it is not unlikely 
that such experiments could be realistically performed in 
the not too distant future, given the rapid rate of im¬ 
provement of isolation from their environment [39[ . 

In general, other properties of the system can also be 
inferred by making preparations and measurements in 
different bases. We list the full set of obtainable quanti¬ 
ties in Appendix D, where we show that when the system 
is highly disordered, the whole spectrum of the operators 


Aq and Ai can be identified. When there are no degen¬ 
erate energy splittings, all diagonal elements of pg in the 
basis of both Aq and Ai can also be determined. 


Vibrational modes 

Another possible application of the formalism pre¬ 
sented here can be found in the context of molecular junc¬ 
tions, where electrons tunnelling between two nano-scale 
electronic binding sites (which could be small molecules 
or quantum dots) couple to the vibrational degrees of 
freedom of a molecule placed between them [43[. The 
binding sites can be modelled as an electron donor and 
acceptor (da) pair - a two level quantum system whose 
energy splitting can be tuned by varying the potential 
of leads to which they are connected; the intervening 
molecule then plays the role of the unknown system S. 
Following Ref. [5[ , we use a spin-boson model to describe 
the da pair interacting with a molecule that has, along 
with its wider environment, a (continuous frequency) col¬ 
lection of vibrational modes: 

H =l\i\d){d\-\a){a\) + V{\d){a\ + \a){d\) 

+ '^ildkld^idl + 7afc|a)(a|) (g) {bk + bl) 

k 

+ <^kb\bki (19) 

k 

where the typical energy gap A (~ 100 — 200 meV) is 
much larger than the tunneling energy E (~1 meV) be¬ 
tween donor and acceptor. To diagonalize the spin part 
of the Hamiltonian we use a polaron transformation [44[ 
generated by the operators' = J2k('‘^dk\d){d\—Uak\0'){a\)'^ 
{bl - bk), where uxk = lxkl<^k- 
After performing the transformation, the energy gap is 
decreased by an amount ~ '“afe)) known 

as the reorganisation energy, and the probe couples to the 
system through a ax operator. If we additionally assume 
that the system and reservoir are both in a thermal state, 
we can then use Eq. (7) to write a simple formula for 
the transition probability from donor to acceptor (see 
Appendix E): 

E 2 

Va-.d [2 - cos (Ar - Ext) exp{/(T)}], (20) 

where 


f{T) = -^(Udj - Uaj)^ COth (1 - COSUJjT), 

in which T is the temperature. 

By measuring paui at different times and for different 
values of A, we can both recover the reorganisation en¬ 
ergy Er - a measure of the overall coupling strength be¬ 
tween V and S ~ and reconstruct the function /(r) in 





Eq. (20). Fourier transforming this function reveals in¬ 
formation about the couplings {'^kd ~ 7fea)^> which for 
realistic systems can be written as a continuous distri¬ 
bution = J2khkd - 7fca)^<5(w - Wfc), known as the 
spectral density. Explicitly, 

/(„) = £Meo.h(^), (21) 

for positive lo. Not only can prominent vibrational modes 
of an unknown system S be identified by reconstructing 
but for a system with known spectral density, /(w) 
can be used to determine the temperature T. That is, in 
this case, our scheme can be used to do thermometry. 

In practice, the range of times for which one could re¬ 
construct /(r) is limited by the constraint in Eq. (13). 
For the parameters mentioned above, the maximum evo¬ 
lution time is ~ 300 ps - much longer than the vibra¬ 
tional coherence time for some typical spectral densities 
in large biomolecules [45]. If necessary, this time could be 
extended by increasing A, since r^ax ~ however, 

this would reduce the transition probability Pa-.dt which 
would then require more experiments to accurately mea¬ 
sure. It should also be noted that other environments 
may contribute to the probe signal even in the absence 
of a sample system. For example, the leads which con¬ 
trol the bias across the da-pair can be considered as a 
fermionic bath [46]. 

V. DISCUSSION 

We have presented a general scheme for using a two- 
level system to probe an unknown system to which it is 
coupled: the application of a precisely controlled external 
bias to the probe, in conjunction with measuring transi¬ 
tion probabilities, allows for the reconstruction of many 
properties of the system, including information about its 
state and spectrum. Importantly, much of this informa¬ 
tion can be obtained with the probe being prepared and 
measured in a single preferred basis {[tto) , [tti)}. 

Since our scheme rests on only a few assumptions - 
notably we do not require the coupling strength to be in 
any particular regime, so long as the probe splitting can 
be made sufficiently large - we expect it to find applica¬ 
tion in many physical scenarios. Our method can be seen 
as a kind of generalised spectroscopy, in that it gives us 
access to the spectrum of an unknown system. In par¬ 
ticular, it could be used to determine the electronic and 
vibrational spectra of large molecular complexes, an im¬ 
portant step in understanding many biological processes, 
such as photosynthesis [47] . 

Unlike traditional spectroscopic techniques, however, 
the probe coupling in our approach is strictly off- 
resonant. In addition to satisfying the conditions in 
Eqs. (12) & (13), which ensure the probe-system inter¬ 
action acts perturbatively, the probe’s splitting must be 
larger than energy gaps between system states which it 
directly couples. 


Though the degree of control required to observe the 
oscillations in probability, and hence probe the system, 
is high, it is certainly within the realms of experimental 
possibility for some systems. In particular, ultrafast spec¬ 
troscopy routinely investigates evolution on timescales of 
order ~ ns with a resolution of ^ fs or even ~ as [48[. A 
comparatively long evolution time of Ins requires control 
over the applied bias on the order of ^ /xeV - though, 
depending on the probe-system coupling, the magnitude 
of the bias may need to be much larger (see Appendix B) . 
For a magnetic field splitting a hyperfine level of an al¬ 
kali atom, this corresponds to a field strength of order 
~ lOOmT. Longer evolution times require an even greater 
degree of control, and larger values of A - the maximum 
value of which may be fundamentally or technically lim¬ 
ited in some systems - to observe oscillations, but for 
many interesting physical processes, coherent quantum 
dynamics occurs on sufficiently short timescales. 

An interesting potential extension of our results would 
be to use the qubit probe to control a more complex sys¬ 
tem S. For example, by regularly removing excitations 
from the probe (using some external time-dependent con¬ 
trol), the latter could act as an energy sink for the system. 
While the interaction with the probe would allow it to 
escape ‘dark’ subspaces (from which it cannot be directly 
cooled), the small effective coupling in the large-A regime 
would prevent it from being directly affected by the ex¬ 
ternal control on the probe. The net effect would be to 
cool the degrees of freedom belonging to S. 

Finally, there are analogies but also important differ¬ 
ences between the ideas outlined in this article and quan¬ 
tum state and process tomography. In quantum state to¬ 
mography [49] , one endeavours to make statements about 
the state of a given quantum system by performing (gen¬ 
eralized) measurements on the latter. These measure¬ 
ments are done by a classical apparatus or a field. On the 
other hand, in quantum process tomography [9, 50, 51] 
known quantum states are used to probe a quantum pro¬ 
cess and describe the quantum dynamics. In the ap¬ 
proach outlined in this article, we are not varying the 
input and the output states, as one would do in a to¬ 
mographic setting. The initial and final states are fixed, 
as is the probe Hamiltonian. We use the quantum sys¬ 
tem V to probe its environment S and then infer some 
features of the latter from the induced dynamics of the 
former. Interestingly, different features are unearthed at 
different orders in the (inverse) energy splitting of V. 

Clearly, one cannot expect a finite quantum probe to 
be able to unveil the full complexity of a much larger 
system with which it interacts. Nevertheless, memory 
effects, energy structure and coherent features do come 
to light, being associated with different couplings and 
different choices of initial states. It would be interesting 
to understand whether an exact solution of the problem 
(say, at all orders in the bias A, if perturbation theory 
is consistent) would be able to probe different (deeper) 
characteristics of a complex system, including its most 
elusive quantum aspects. 
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Appendix A: Change of basis of control 

The system operators for different choices of the control parameters can be related with reference to their counter¬ 
parts in the z-direction (which defines 0 and (j)) as 

Aq =cos^ 9 Aq + sin^ 0A^ + cos0sm0(B^ + (Al) 

A‘[ = sin^ 0Ag + cos^ 0A^ - cos 0 sin 0(B^ + B^^), (A2) 

=e-^^[sm^0B^'^-cos^0B^+cos0sm0(A§-Af)]. (A3) 


Appendix B: Conditions for perturbation series to exist 


In order to only consider contributions at first and second order in 1 /A, we require all higher order terms in Eq. (7) 
to be strictly smaller in magnitude. Here we specify in detail the conditions on A for this to be the case. 

Each term in Eq. (6) can be written in the form 


(zA)^ 


where = f dsi f ds 2 ■ ■ ■ f ds, 

Jo Jo Jo 


X ^Si 77^2 


■H, 


(Bl) 


and r = x + y. Taking the trace over these terms, with the relevant projection operators for preparation and 
measurement states, leads to the following contributions to pp-,ci'. 


(zA) 


-tr 


1 O 


^ ^ a:a^b;b, (j„| PS |fc™) l,) l^\ C/(-)|7r„ Kp Q - (^2) 


(zA) 


nmpq 

jkll' 


where we have expanded the system operators in their eigenbases as Aq = J2j ^j\Jo){jo\ and Ai = ^j\ji){Ji\- 

Therefore, a valid perturbation expansion requires A S> | (tt, lq\ Ut^'^ {iTn jn){T^m \tTp Ip) ^ can now 

use the definition of in Eq. (Bl), as well as the interaction Hamiltonian given in the main text to write these 
amplitudes explicitly. Firstly, note that 


TTqlq\ \TTn jn) OC 


^n,q 

1 


X even, 
X odd. 


(B3) 


We will therefore treat these two cases separately. For the even case, we have 




fc(l)fc(2)...fe(^-l) 
(‘Xt 
X 


(a:-l)\/ Ax-l) 




Jo 


[ ["ds,-- - 

Jo Jo Jo 

(l+(S“(_2 ) )/a) g—. + )/a) 


(B4) 
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and 


(^1 h 1ki Ji) = 5] {h B fc W )(fc W i?t fcf))(*(2) 


fc(l)fc(2)...fe(^-l) 
[‘Xt 
X 


(a;-l)\/ 7 (x-l) 




Jl 


[ "dsi Tds^... 

^0 -/o -/o 

^ ^ ^ (l+(£”(._i) )/a) (l + (i5:(.-i) -k)/^) . 


(B5) 


For the odd case: 


(^1 h 1ko jo) = Y1 (^1 ^ 


fc(l)fc(2)...fe(--l) 
[‘Xt 
X 


(a;-l)\/ 7 (a:-l) 




Jo 


[ "dsi Tds^ ... 

^0 -/o -/o 

^ ^ ^ ZS_1 (i+(^J:(.-2, )/a) (l + (i5«-i5^(_i, )/a) ^ 


(B6) 


We can use these expressions to relate the magnitude of lambda required for a good perturbation series to the operators 
Ak and B. In order for the contribution from the term proportional to 1/A’’ to be smaller in magnitude than that 
proportional to l/A*", if we ignore the time-dependent part it is sufficient that 


> I (jil B ko) I, 


(B7) 


for all |ji) and ko) reachable from ps within w -|- 1 applications of the interaction Hamiltonian (i.e., 
(ji(A:o)| Iji(fco)) k 0), and 


A>|(ji|.Sko) 


(B8) 


for all l/i) and ko) reachable from pg within r interactions. 

When A does not equal any of the energy splittings in the system (we will consider the effect of resonances below), 
iteratively performing the first few nested integrals leads to 

X'dsi rds2 ... 

Jo Jo Jo 

Jo 


dsi 


ds2 ... 


^zs_i {(eI^^-2) -£;)(.-!) )-(£;“-£;)(._i) ))/a _ zs.-i (i-hs/(.-2, )/a) 


pXr 




dsi 


ds2 ... 


1 + {E^-El^_,0/X 

Jo 


X A- 


-£;)(.-3) )/A-(i3“(^_2) )/A-H£;“-iS)(._i, )/a) _ -is,-2 -£^((.-3) )/a) 


m,^-2, - - (E^ - El^_,0) (1 + (E^ - El.-^O/^) 

^-Z3„_2 ((£“(_„/A _ ^-Z3„_2 (!+(£“(_,,-£;)(_3, )/a) 


(1 + (ifO - El^_, 0 /X) (1 + (i?,V-2, - El^_, 0 /X) 


(B9) 


As can be seen in the second to last line of Eq. (B9), factors of A appear in the numerator whenever the I’s cancel in 
the exponent (one can see that factors of 1/A in the denominator contribute to higher order terms in the expansion). 
For a given term, this will happen at most every other integration (e.g. when integrating over Sa;_i, Sx- 3 , Sx -5 etc.), 

X . X —1 

hence the leading order term of the whole integral will be explicitly proportional to A 2 in the even case and A “ 2 “ in 
the odd case. If we were to require the magnitude of such terms to still be strictly smaller than those at lower orders. 
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we would need that the ratio of subsequent terms goes as 1/A. Ignoring the complex exponential factors resulting 
from the integration, this constraint can be written as 


-V£, E,, and Vx < n/2 (n even), 

|tr{(EBt).|i„)(i„|gBt||;)(,;|(BBt)(e«)/2-e,,^)|^|Ee _ e.| ^ 1 „ „ 

~A^ t., and Vx < r/2 (r odd), 


(BIO) 


and similar expressions with additional pairs of projectors sandwiched between B operators and corresponding factors 
of l/|if° — E^i\- The same should also hold with 0 —>■ 1 and B ^ BB^. 

There will also be terms appearing in Eq. (B9) where the pairs of energy levels in the exponent cancel, i.e., 
{El — El) — {El, — El) = El — E^, = 0. Even when there are no degeneracies, this can occur when I = V. These 
terms introduce extra factors of At, leading to a restriction on the time for which the perturbation series remains 
valid. Interacting for a time r places constraints on A identical to those in Eq. (BIO), but with l/|i?/ — Ei,\ replaced 
by T (and energies required to be equal, Ef = Ei,). 

However, for the results in the main text to hold we do not need each subsequent term in Eq. (6) to be strictly smaller 
than the last, only that later contributions are much smaller than the zeroth, first or second order terms we present 
in Eqs. (8-11) (depending on the choice of preparation and measurement). For the zeroth order case, in addition 
to Eqs. (B7) & (B8), we must also worry about contributions from second order terms of the sort discussed below 
Eq. (B9). Though they are nominally proportional to 1/A after taking into account the integration, the coefficients 
could in principle be large; considering expressions similar to those in Eq. (BIO), where the second order terms are 
required to be much smaller than those at zeroth order, it is straightforward to show that sufficient constraints on A 
are 


A»5]4,|(Zo|Bti?|Q(^(,|P5lMI where 4 , = 

w El = Ef, 

and the same with 0 —)■ 1 and B^ B ^ BBK 

For the first order terms in the main text to accurately describe the transition probabilities, we require that the 
second order terms are much smaller. This leads to A-independent limits on ps and the Ak and B operators: 

tr{Hp5} » B^B\1'Mps |/o) I, (B12) 

w 


and the same with 0 —)■ 1 and B^B ^ BB\ where k/j, is defined as in Eq. (Bll). 

When both preparation and measurement are in the control basis, the corresponding constraints are again dependent 
on A, since for the case ()3| a) = 0 (trivially related to the case with (/3| a) = 1), the next highest terms are at fourth 
order in the perturbation series (r = 4 in Eq. (B2)). Strictly they are 


, V- 0 I ('ol |io> 

-|,r{B.pxB)|- 


(B13) 


and the same with 0 —1 and B ^ B^. 

In the case where there are energy differences Ej — El = A, then, in general, extra powers of At appear 
from the nested integrals, and the perturbation expansion is no longer valid. However, if we restrict such en¬ 
ergy differences to transitions which require at least w interactions to reach from the initial system state, i.e., 

I (/o| {HsYps{Hs)^ |jo) I — I (^i| {iis)^ps{Hs)^ \ki) I ~ 0 for z < in, then terms in Eq. (B 2 ) with explicit 1 /A’' de¬ 
pendence will include at most r — 2w additional powers of A. Formally, we require A 7 ^ \Ej — El\ for all |A:i), |/o)j 

such that tT[\jo){jo\Hs\ki){ki\Hg psH^] 7 ^ 0 , where x + y < w — 1. 

For continuous frequency systems (often the case when a system of interest cannot be isolated from its wider 
environment), such resonances are unavoidable; A = Ej — El always, for some E^j and El- The couplings for such a 
system can be rewritten in terms of continuous functions G{uj) = /c (^il l/o) ~ \Ej — El\). Our perturbation 
expansion still holds when G(A) < 1/(A^t) - effectively removing the offending terms to leading order. 

With general |a) and |/3), all terms of the form given in Eq. (B2) contribute to pp-.a, so the leading order contribution 
from the 1/A’' term is actually A’’’/^A^/^/A’' = 1/A’'/^ (x and y both even), A*^’’’“^^/^A*^/^/A’' = l/A^’'"*'^^/^ (one even, one 
odd) or aO-i)/2aO-i)/2/;s^’' = l/A’'/ 2 +i (x and y both odd), depending on the value of r (and neglecting the resonance 
conditions described above). After those terms treated explicitly in the main text, terms with r = 1 and r = 2 have 
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the largest contributions, but neither has contributions at leading (zeroth) order, thus they can be neglected. This 
remains the case as long as there are no energy gaps resonant with A within w = 1 interaction of ps. 

In the case where |a) or |/3) is in the control basis, we must consider whether terms with r = 2 or greater contribute 
at first order in 1/A. Using Eq. (B9), we find that r = 2 and r = 3 terms are both proportional to 1/A^ at leading 
order, as long as there are no w = 1 resonances. 

Finally, when |a) and |/3) are both in the control basis, we need to make sure r > 3 terms do not contribute at 
second order. Though r = 3 terms go to zero in this case, we might worry that for r = 4, terms with a; = 2, j/ = 2 in 
Eq. (B2) do have second order contributions. However, when (a| /?) = 0 (which is trivially related to the case where 
they are identical), Eq. (B3) implies that such terms disappear; therefore, without resonances, r = 4 contributes at 
third order at most (with a; = 1 and y = 3) and r = 6 at fourth order. In order for there to be no contributing 
resonances, we require that A is greater than all energy differences within w = 3 interactions. 

To summarise: for the results in the main text to hold, we firstly require that A (jo| I? |fci), at least for those 
energy levels |jo) and \ki) which are within three interactions of the initial state ps {{jo\ BB^BpgBB^B \ki) ^ 0). 
We also require that all other matrix elements (jo| B \ ki) are not significantly greater than A. Furthermore, we need 
all energy differences that involve aforementioned energy levels, reachable within three interactions, to be different 
from the probe splitting: A ^ \Ej — E/|. There is a final requirement given by Eqs. (B11~B13), depending on the 
leading order term of pp-a (and hence choice of preparation and measurement) that relates A to the interaction time 
r and the properties of the system for zeroth or second order leading 1/A dependence. For first order dependence, 
there is a fundamental constraint on the properties of the system. 


Appendix C: A-dependence of 


Order one — By expanding the system operators in their eigenbases as Ag = and Ai = 

Sn *^an perform the integrals in Eq. (6). For the first order term, we have 


ds |tr[aoai (cr(e_0))^e*®H„P5 - &o6ie ^^Hr[(|aipp5i?J - |aopHl/ 05 )e*^“^e | + c.c. 

1 1 _ / . n 1 

= - (toi|H^ |no) (?^o| ^i) {h\ Ps |mi) 


A e ^ + iE^n- El)/X 

|2 /„ I I, \ /, I _ \ I Dt 


- hlhi ^ |ao|^ (nol ps Xo) (^o| wi) (mi| B^ |no) - oga* {cr(e, 4 ,))p {no\B\mi){mi\ps |ng) j + c.c. 

i X 

y {aoal (^(S.0))^ ( (S)5 - e*^") + hlh^\a^\^ B^^ B’^)^ e*^") 

+C.C. + 0 , (Cl) 


where in the final line we see and as defined in Eqs. (8) & (9) appearing. The second order 

(and higher) terms that result from expanding the denominator in the second line are absorbed into the definition of 
in Eq. (7). 

Order two — Similarly, we can perform the integrals in the second order terms in Eq. (6). We ignore those terms 
which never appear at leading order in 1/A, which are included in the definition of /'^^^(A). For the remaining terms 
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System quantity Fourier transform 


1 

2 

3 

4 

5 

6 

7 

8 
9 




(nil p5 |mo) (mol ni) (5(a; - {Ei - E^)) 
(no| B|mi)(mi|p 5 |no) 5(w - {E° - E^)) 
(mil no) (noI 5^5 |mi) 5(a; - {E° - Ela)) 
(no| mi) (mil psB"' |no) 5(a; - {E'i - El^)) 
(nfel B |mfc) (mfel ps \nk}S{oj - {E^ - E'^)) 


E™ K| B|mi)(mi|Btp5 |no) 3 {lo - (F° - F^)) 
(i3te*^«"Be-*-"iEs E„^ («o| B|mi)(mi|p5B^ |no) 5(u; - {E° - E^)) 
(e‘-^«"BBte-‘^«Es E„™ (^o| BBVo)(mo|pi5 |no) S{uj - {E° - El)) 
(e--^^"i?tj3e-^^E5 E.^ (nil i3tB|mi)(mi|pi; |ni) - [Ej - El)) 


Significance 


Equal to 1 when j4o = Ai = Hs. 


Witness for initial system coherence. 
Two-time correlation functions. 


TABLE II. Properties of the system deducible from probability oscillations in A. Fourier distributions also presented. 
Here we have written the environment operators in their eigenbases as = E„ |wfe)(nfe|. 


we find 


pAr pAr 


/ ds I ds'tr 


/o Jo 

" A2 2 


^ f |ao|2 (nol i3bi)(pi|BVo)(mo|p5 |no) ds ds' 


- |aip (nil p5|mi)(mi|B'l’|po)(Pol-B |tii) [ dse*^^ ^p)/^)® f ds'e K)/^> 

Jo Jo 

- Imp {{B^B)^ + ^ + O , 


A2 2 


(C2) 


where again, the denominator has been expanded in powers of 1/A after integrating. In the last line, we see the terms 

/Q\ / 2^ 

corresponding to '(A) and '(A) defined in Eqs. (10) & (11) respectively. 


Appendix D: Summary of obtainable properties of the system 

Table II shows a summary of all the information one can obtain about the environment by measuring the amplitude, 
phase and displacement of probability oscillations in A as a function of r. The functions have been Fourier transformed 
to give a distribution in energy; for a system with a discrete spectrum, each peak in the distribution will correspond 
to pairs of eigenstates of Aq or Ai. In general, there will be multiple such pairs contributing to each peak, but for 
highly disordered environments (randomly distributed energy levels), each pair of levels will have its own peak. 

When no two energy differences in the system are the same, one effectively has access to each of the terms in the 
energy distributions listed in Table II. Furthermore, the whole spectrum of the operators Aq and Ai can be identified 
by comparing the different frequencies present in the distributions in row 5; triplets of energy levels can be identified 
from pairs of frequencies which sum to a third (e.g. (£!„ — E^n) + {Em — Ep) = E^ — Ep). By finding all such 
triplets which include the lowest and highest energy states (corresponding to the highest frequency present in the 
distribution), the entire spectrum can be reconstructed. 

With information of the spectrum, one can pick out all the terms in, for example, the first row which come from 
the same energy level E^ The sum of these terms gives the diagonal elements of pe in the basis of Aq, (no| ps jn-o). 
Similarly, the diagonal elements in the basis of Ai can also be found. Weighting the diagonal matrix elements by 
the corresponding energy eigenvalues allows the expectation values (Aq) and (Ai) to be calculated. Importantly, the 
sum of these gives the mean energy of the environment (He) = (Aq) + (Ai). However, it should be noted that this is 
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only possible when there is coherence in the initial system state - something that it is always possible to achieve by 
rotating the control basis, i.e., choosing different values of 9 and cj). 


Appendix E: Vibrational modes: derivation of probabilities 

Here we outline the details of the derivation of the probabilities in subsection IV. We start with repeating the 
Hamiltonian of the molecular junction spin-boson model: 

H =^A (|d)(d| - |a)(a|) + V{\d){a\ + |a)(d|) -f '^{'ydk\d){d\ -f 7afc|a)(a|) (g) {bk + bl) + '^uJkblbk, (El) 

k k 

with the relevant parameters explained in the core text - note that, for convenience, we have set the zero of energy 
to be halfway between the donor and acceptor states of the probe. Next, we apply the polaron transformation 
H —>• e^He~^ , generated by: 


s = '^{udk\d){d\ - Uak\a){a\) 0 {bl - bk): (E2) 

k 

This transformation leaves operators \X){X\^^^ ^ unchanged but modifies the Hamiltonian (up to a constant energy 
shift) to: 

H Er){\d){d\-\a){a\) + H-ps, with (E3) 

Hps = M)(a| ® B + |a)((i| ® + \d){d\ (g) Hq -|- |o)(a| (g) Ai, (E4) 

where Wk ~ '^ak) written explicitly, the operators Hop and B relevant for our formalism are 

Ho = Hi = ^ ujkblbk, (E5) 

k 

B = VIlkD{udk - Uak)- (E6) 


In the above, the displacement operator is defined as D{^) = exp{^6^ — ^*b}. A natural basis for preparations and 
measurements is the basis \X){X\^^^ As a consequence, only the second-order terms play a role in the analysis of 
the probabilities pp-.a- Since we have B^ B = 1, the probabilities simplify to: 


pp-.c =qp-.o. - - |aip)Re {2V^ - {B ^, (E7) 

where qp:a G {0,1}, (|aoP — |aiP) = ±1 and {az)p = ±1, depending on the values of a and /3. Restricting ourselves 
to transitions from d to a and substituting in the definition of B, we can further write: 


Pa:d X [2 - . 


(E8) 


To analyze the constraints on the applicability of the method in this scenario, we first consider a transition element 
|(mim 2 • • • TOfc • • • |i?|nin 2 • • • rife • • • )|, where the fcth mode initially has occupation number rife and ends up with 
occupation number nife. This quantity must be much smaller than A for states with non-negligible initial probability 
(and not too large for other states). Using the well known expression for the matrix elements of displacement 
operators [-52], we can write the constraint on A as 


A » Ve-5 —)^nfe (' \udk - Uafer-^A™-” {Wdk - UafeH , (E9) 

where L^{x) is a Laguerre polynomial. Unless there are particularly strongly coupled low frequency modes (corre¬ 
sponding to a sub-Ohmic spectral density, for which the polaron transformation is known to have problems), the 
satisfaction of this inequality principally depends on the value of V, which we have chosen to be small. The second 
condition related to the validity of our expansion is given in Eq. (B13). Taking advantage of the fact that in this 
example B^ — B~^ we can easily show that the relation between A and r simplifies to: A rU^. 
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With the additional assumption on the form of system and the reservoir states being both thermal, we can evaluate 
the expectation value of the displacement operators to get the suppressing factor: 

(UjD {{udj - -!))> = exp{-Ej(^^<t “ - cosuJjt) coth 

with Pb being the inverse temperature Boltzmann factor. This allows us to rewrite the probabilities ppia in a much 
simpler form presented in the main text in Eq. (20). 
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